
#include <fcntl.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <CL/cl.h>
#include <jni.h> 

const char *KernelSource = "\n" \
"__kernel void square(                                                   \n" \
"   const unsigned int slice,                                  			 \n" \
"   const float omega,                     				             	 \n" \
"   const unsigned int array_rows,     		                             \n" \
"   const unsigned int nthreads,                                       	 \n" \
"   __global float* y                                                    \n" \
")                                           						   	 \n" \
"{                                                                     	 \n" \

"int id = get_global_id(0);												 \n" \
"if(id<nthreads){														 \n" \
"	int ilow = id * slice;												 \n" \
"	int iupper = (id+1)*slice; 											 \n" \
"	if (iupper > array_rows ) iupper=array_rows; 						 \n" \
"	int i;																 \n" \
"	float x;              												 \n" \
"	float x0 = 0.0;       												 \n" \
"	float x1 = 2.0;        												 \n" \
"	int nsteps = 1000;     												 \n" \
"   float dx;															 \n" \
"   float rvalue;        												 \n" \
"   int select = 2;        												 \n" \
"   int j;	        													 \n" \
"	if(id==0) {															 \n" \
"		ilow=id*slice+1;												 \n" \
"		x = x0;															 \n" \
"		dx = (x1 - x) / nsteps;											 \n" \
"		rvalue = pow(x0+(float)1.0,x0) / (float)2.0;					 \n" \
"		if (nsteps != 1)												 \n" \
"  		{																 \n" \
"			--nsteps;               									 \n" \
"			while (--nsteps > 0)										 \n" \
"			{															 \n" \
"		    	     x += dx;												 \n" \
"                            rvalue += pow(x+(float)1.0,x);							 \n" \
"                       }															 \n" \
"		}																 \n" \
"		rvalue=(rvalue + pow(x1+(float)1.0,x1) / (float)2.0) * dx;	     \n" \
"   	        y[0] =(rvalue/ (float)2.0);										 \n" \
"	} 																	 \n" \
"	for (i=ilow; i<iupper; i++) {										 \n" \
"		x = x0;															 \n" \
"		nsteps = 1000;													 \n" \
"		dx = (x1 - x) / nsteps;											 \n" \
"		rvalue=(pow(x0+(float)1.0,x0)*cos((omega*i)*x0))/(float)2.0;	 \n" \
"		if (nsteps != 1)												 \n" \
"  		{																 \n" \
"			--nsteps;               									 \n" \
"			while (--nsteps > 0)										 \n" \
"			{															 \n" \
"		    	    x += dx;												 \n" \
"                           rvalue +=(pow(x+(float)1.0,x)*cos((omega*i)*x)); 	 	 \n" \
"                       }															 \n" \
"		}																 \n" \
"		rvalue=(rvalue + (pow(x1+(float)1.0,x1)*cos((omega*i)*x1))/(float)2.0) * dx; \n" \
"   	        y[i] = rvalue;													 \n" \
"		x = x0;															 \n" \
"		nsteps = 1000;													 \n" \
"		dx = (x1 - x) / nsteps;											 \n" \
"		rvalue=(pow(x0+(float)1.0,x0)*sin((omega*i)*x0))/(float)2.0;	 \n" \
"		if (nsteps != 1)												 \n" \
"  		{																 \n" \
"			--nsteps;               									 \n" \
"			while (--nsteps > 0)										 \n" \
"			{															 \n" \
"		    	    x += dx;												 \n" \
"                           rvalue +=(pow(x+(float)1.0,x)*sin((omega*i)*x)); 		 \n" \
"                       }															 \n" \
"		}																 \n" \
"		rvalue=(rvalue + (pow(x1+(float)1.0,x1)*sin((omega*i)*x1))/(float)2.0) * dx; \n" \
"   	        y[array_rows+i] = rvalue;										 \n" \
"	}						 											 \n" \
"}                                                                       \n" \
"}                                                                       \n" \
"\n";



////////////////////////////////////////////////////////////////////////////////
JNIEXPORT jfloatArray JNICALL Java_section2_series_SeriesOCLTest_seriesFromC
  (JNIEnv *env, jobject jobj, jint jslice, jfloat jomega, jint jarray_rows, jint jnthreads)
{
    int err;                            // error code returned from api calls
     
    jfloatArray res;          			// results returned from device
   
    size_t global;                      // global domain size for our calculation
    size_t local;                       // local domain size for our calculation
 
    cl_device_id device_id;             // compute device id 
    cl_platform_id platform;		// compute platform id
    cl_context context;                 // compute context
    cl_command_queue commands;          // compute command queue
    cl_program program;                 // compute program
    cl_kernel kernel;                   // compute kernel

	cl_mem y;

	unsigned int nthreads = (unsigned int)jnthreads;
	unsigned int array_rows = (unsigned int)jarray_rows;
	unsigned int slice = (unsigned int)jslice;
	float omega = (float)jomega;
	res = (*env)->NewFloatArray(env, 2*array_rows);
    float* result = malloc(2*array_rows * sizeof(float));
    
    // Connect to a compute device
    //
    int gpu = 1;
    err = clGetPlatformIDs(1, &platform, NULL);
    if (err != CL_SUCCESS)
    {
        printf("Error: Failed to select paltform!\n");
        return res;
    }

    err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device_id, NULL);
    if (err != CL_SUCCESS)
    {
        printf("Error: Failed to select device!\n");
        return res;
    }

    //
    context = clCreateContext(0, 1, &device_id, NULL, NULL, &err);
    if (!context)
    {
        printf("Error: Failed to create a compute context!\n");
        return res;
    }

    // Create a command queue
    //
    commands = clCreateCommandQueue(context, device_id, 0, &err);
    if (!commands)
    {
        printf("Error: Failed to create a command queue!\n");
        return res;
    }

    // Create the compute program from the source buffer
    //
    program = clCreateProgramWithSource(context, 1, (const char **) & KernelSource, NULL, &err);
    if (!program)
    {
        printf("Error: Failed to create compute program!\n");
        return res;
    }
 
    // Build the program executable
    //
    err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
    if (err != CL_SUCCESS)
    {
        size_t len;
        char buffer[2048];
 
        printf("Error: Failed to build program executable!\n");
        clGetProgramBuildInfo(program, device_id, CL_PROGRAM_BUILD_LOG, sizeof(buffer), buffer, &len);
        printf("%s\n", buffer);
        exit(1);
    }

    // Create the compute kernel in the program we wish to run
    //
    kernel = clCreateKernel(program, "square", &err);
    if (!kernel || err != CL_SUCCESS)
    {
        printf("Error: Failed to create compute kernel!\n");
        exit(1);
    }

    // Create arrays in device memory for our calculations
    //
    y = clCreateBuffer(context,  CL_MEM_WRITE_ONLY,  sizeof(float) * (2*array_rows), NULL, NULL);

    if (!y)
    {
        printf("Error: Failed to allocate device memory!\n");
        exit(1);
    }    

    // Set the arguments to our compute kernel
    //
    err = clSetKernelArg(kernel, 0, sizeof(unsigned int), &slice);
    err |= clSetKernelArg(kernel, 1, sizeof(float), &omega);
    err |= clSetKernelArg(kernel, 2, sizeof(unsigned int), &array_rows);
    err |= clSetKernelArg(kernel, 3, sizeof(unsigned int), &nthreads);
    err |= clSetKernelArg(kernel, 4, sizeof(cl_mem), &y);
    if (err != CL_SUCCESS)
    {
        printf("Error: Failed to set kernel arguments! %d\n", err);
        exit(1);
    }
 
//get the number of compute units (SMs) of this device
    cl_uint units;
    clGetDeviceInfo( device_id, CL_DEVICE_MAX_COMPUTE_UNITS,  sizeof(cl_uint),  &units, NULL);

//get the multiple of local work group
    size_t mult;
    err = clGetKernelWorkGroupInfo(kernel, device_id, CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE, sizeof(mult), &mult, NULL);
    if (err != CL_SUCCESS)
    {
        printf("Error: Failed to retrieve kernel work group info! %d\n", err);
        exit(1);
    }
  
    global = nthreads;
    local = global / (units); // create at least 1 warp per multiprocessor (a group has at least one warp)
    if(local>mult) {
    local = local / mult; //now make local a multiple of mult (warp size) so that resources in the warp are not wasted
    local = local * mult;
    }

    
    if(global % local != 0)
    {
		global=global-(global%local);
		global+=local; //make global a multiple of local
    }

 
 	struct timeval tim;
	double starttime, endtime,temp;
	//measure kernel execution time
   	gettimeofday(&tim, NULL); 
	starttime =  tim.tv_sec*1000+(tim.tv_usec*0.001);

    err = clEnqueueNDRangeKernel(commands, kernel, 1, NULL, &global, NULL, 0, NULL, NULL);
    if (err)
    {
        printf("Error: Failed to execute kernel!\n");
        return res;
    }
    
        
	
    // Wait for the commands to get serviced before reading back results
    //
    clFinish(commands);
    
    gettimeofday(&tim, NULL); 
	endtime  = tim.tv_sec*1000+(tim.tv_usec*0.001);
	printf("%4.3f", (endtime-starttime)*0.001);
    // Read back the res from the device
    //
    err = clEnqueueReadBuffer( commands, y, CL_TRUE, 0, sizeof(float) * (2*array_rows), result, 0, NULL, NULL );  
    
    if (err != CL_SUCCESS)
    {
        printf("Error: Failed to read output array! %d\n", err);
        exit(1);
    }
    
    // Shutdown and cleanup
    //
    (*env)->SetFloatArrayRegion(env,res, 0, 2*array_rows, result);


	free(result);
    clReleaseMemObject(y);
    clReleaseProgram(program);
    clReleaseKernel(kernel);
    clReleaseCommandQueue(commands);
    clReleaseContext(context);

    return res;
}
